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ABSTRACT 

A  general  uniform  asymptotic  expansion  is  constructed  to  remove  certain 
anomalies  which  appear  in  the  ray  method  expansion  for  unsteady  flow  in 
unconfined  aquifers,  and  makes  the  asymptotic  method  based  upon  the  principles 
of  geometrical  optics  even  more  useful.  In  particular,  expansions  uniform  in 
a  region  containing  a  caustic,  a  line  of  zero  depth,  a  combination  of  both, 
two  caustics,  and  two  lines  of  zero  depth  are  explicitly  constructed. 
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UNIFORM  ASYMPTOTIC  EXPANSIONS  FOR  UNSTEADY  FLOW 
IN  UNCONFINED  AQUIFERS 

M.  C.  Shen 

I.  Introduction 

The  problem  of  fluid  flow  in  unconfined  aquifers  has  attracted  much 
attention  in  recent  years,  because  of  its  great  importance  in  underground 
water  flow,  land  subsidence,  oil  flow  to  wells,  and  many  others.  Several 
approximate  theories  for  the  study  of  such  a  problem  are  now  at  our  disposal, 
for  example,  the  linear  equations,  the  Boussinesq  equation  (Scheidegger  [1]), 
and  the  delayed  field  equation  (Boulton  [2],  Herrera  et  al  [3]).  However, 
these  equations  are  generally  used  to  deal  with  homogeneous  aquifers  of 
uniform  depth  only.  In  a  previous  paper  (Shen  [4]  ),  we  developed  an 
asymptotic  method  employing  rays  as  in  geometrical  optics  for  the  study  of 
underground  flow  in  unconfined  nonhomogeneous  aquifers  of  variable  depth 
governed  by  linearized  equations.  It  was  motivated  by  the  work  on  diffusion 
equation  (Cohen  and  Lewis  [5],  Voronka  and  Keller  [6]),  and  on  surface  waves 
on  a  fluid  of  variable  depth  (Keller  [7],  Shen  [8]).  The  asymptotic  expansion 
for  the  solution  of  the  problem  consists  of  a  phase  function  and  successive 
approximations  to  an  amplitude  function.  The  phase  function  satisfies  the 
Hamilton-Jacobi  equation  and  may  be  solved  by  the  method  of  characteristics, 
which  in  turn  determines  a  family  of  curves  called  rays.  The  successive 
approximations  to  the  amplitude  function  satisfy  a  transport  equation  along 
the  rays  and  can  be  determined  by  simple  integration.  However,  the  ray  method 
fails  at  the  places  where  the  rays  form  an  envelope  or  the  depth  of  the 
aquifer  is  zero,  and  the  amplitude  function  becomes  infinite  there. 
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The  purpose  of  this  work  is  to  remove  these  anomalies  by  means  oi  ~iform 
asymptotic  expansions  in  terms  of  solutions  of  general  second  order  ordinary 
differential  equations  in  a  manner  as  suggested  in  Shen  and  Keller  [9] .  By 
appropriately  choosing  the  coefficients  in  the  differential  equation,  and 
using  the  solution  of  the  equation  in  the  asymptotic  expansion,  the  amplitude 
function  then  becomes  finite  at  these  anomalies.  The  refined  method  is  no 
more  difficult  to  use  than  the  ray  method,  since  by  a  simple  transformation 
the  results  in  the  former  are  reduced  to  those  in  the  latter.  The  uniform 
asymptotic  expansion  at  a  caustic  is  a  slight  modification  of  the  method  due 
to  Ludwig  [10].  However,  the  behavior  of  a  solution  at  a  line  of  zero  depth 
of  an  unconfined  aquifer  is  studied  here  for  the  first  time. 

We  formulate  the  problem  and  recapitulate  the  ray  method  in  §2.  In  §3  we 
consider  the  uniform  asymptotic  expansion  under  the  shallow-water  assumption, 
which  is  simpler  to  use.  In  §4,  we  consider  several  special  cases,  for  which 
comparison  equations  are  explicitly  constructed.  Some  of  the  detailed 
derivations  are  given  in  the  appendices. 


2 .  Formulation 


The  linearised  equations  governing  the  compressible  flow  in  an  unconfined 
aquifer  are  assumed  to  be  the  following: 


7*  .  (K*V*$*)  - 

(1) 

S2*t*  +  K*^z*  "  ° 

at  z*  ■  0  , 

(2) 

7*^*  .  7*(z*  +  h*)  *  0 

at  z*  ■  -h*  . 

(3) 

The  physical  meaning  and  derivation  of  (1)  to  (3)  may  be  found  in  Scheidegger 

[1].  Here  is  the  potential,  V*  *  (3/3x*,  3/3y*,  3/3z*),  z*  is  positive 

upward,  K*(x*,y*,z*),  S*(x*,y*,z*) ,  and  S*(x*,y*),  are  given  positive 

functions  assumed  to  be  sufficiently  smooth.  The  free  surface  of  the  flow  is 

z*  “  n#(x*,y*,t*)  “  -♦Vg  at  z*  -  0  where  g*  is  the  constant 

gravitational  acceleration,  and  z*  *  -h*  is  the  rigid  bottom,  nonpositive 

and  sufficiently  smooth.  We  measure  x*,  y*,  z*  and  h*  in  units  of  H, 

the  mean  depth  of  the  aquifer,  t*  in  units  of  (H/g)1^ ,  Assume  that 

K*(x*,y*,z*)  is  bounded  below  by  a  positive  constant  M  and  measure  K*  in 

units  of  M.  Now  let  the  horizontal  length  scale  be  L  and  assume 
3/2 

6  *  L/H  where  0  is  a  large  parameter.  The  nondimensional  variables  are 

x  -  0“3/2xVh,  y  -  0~3/2y*/H,  h  -  h*/H  , 


z  -  z*/H,  K  -  K*/M,  t  -  0-2 tVO/g)^  , 

S|  -  S*H3/2^2/m>  s2  ,  S*(g/H)1/z/M,  ♦  -  ♦V(gH)  . 

In  terms  of  them,  (1)  to  (3)  become 


V  •  (K V*)  +  03(k*  )  -  8$  , 

z  Z  It 

(4) 

S,<t>  +  e2K$  *  0  at  z  -  0  , 

&  t  z 

(5) 

•  Vh  +  82$  ■  0  at  z  *  h  , 

z 

(6) 

where  V  »  {3/3x,  3/3y) . 


Assume  that  (4)  to  (6)  possess  an  asymptotic  expansion  of  the  form 


(7) 


♦  ~  exp[-00 (x,y,t )  ]  (Ajj  +  0_1A1  +  ®”2ft2  +***’  >  • 

Substitution  of  (7)  in  (4)  to  (6)  yields  a  sequence  of  governing  equations  and 

boundary  conditions  for  the  successive  approximations  AQ,  A1#  A~ , . . . .  and 

the  results  are  summarized  as  follows.  Their  derivations  may  be  found  in  Shen 

[4].  Let  ui  *  -0fc,  h  *  and  k2  »  (70)2.  Then  0  satisfies 

2 

u  *  Dk  (8) 


where 

D  -  J®h  K  dz/I,  I  *  S2  +  /®h  S1  dz  .  (9) 

The  characteristic  equations  for  (8)  (Courant  and  Hilbert  [11])  are 


dt/do  =»  p,  dr/do  *  2p  Dk,  du/d  a  »  0 
dk/da  -  -pk27D,  d0/dO  -  pu  , 


(10) 


where  r  =  (x,y)  and  p  is  a  proportional  factor.  The  solutions  of  (10) 
determine  a  two-parameter  family  of  time-space  curves  called  rays,  given  by 

r  -  r [0,ay,a2),  t  -  t{a,a^,a2)  ,  ( 1 1 ) 

where  o^ ,  are  constant  along  each  ray.  By  integrating  d0/do  -  pu  along 
a  ray,  we  have 

0(0)  -  0(o  )  +  J®  pudo  , 

0 

where  6(0^)  is  the  value  of  0  at  some  initial  point  on  the  ray. 

From  the  equations  for  the  second  approximation,  we  obtain 


where 


-S2u  At(x,y,0,t)  -  S2AQt  -  AQK(x,y,-h)k  •  Vh 
-  /°h  UBf  -  Kk2)A1  +  Aq  7  •  Kk  +  2K  k  •  7aq 

A,  -  AQ(-k2L  +  m)  +  f(x,y,t)  , 


siAotIdr 


(12) 


L(x,y,z) 


/!h  K_1(x»y»*>  /*h 


K(x,y,z' )dz*  dz 


I 
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M(x,y,z)  -  K_1(x,y,z)  /*h  S1 (x,y,z* )dz* dz  . 

By  (8)  and  (10),  (12)  can  be  reduced  to  the  following  fora 

(IA^)t  +  V  •  (IAq  dx/dt)  ♦  2NA2  -  0  ,  (13) 

where 

N(x,y,t)  -  -k2[/°h  *  Xk2)L(x,y,z)dz  +  S2«L(x,y,0)] 

+  <*>[ /°h  (S^  -  Kk2  )M(x,y,z)dz  +  S2<»>M(x,y,0 )]  . 

(13)  can  be  integrated  along  a  ray  to  yield 

XA2J(t)exp(-/^  2NI  'dt)  ■  constant  ,  (14) 

where  J(t)  is  the  Jacobian  of  transformation  from  the  ray  coordinates 
o  -  t,  ol#  o2  to  the  coordinates  x,  y,  t.  At  h  ■  0  and  at  a  caustic 
where  J(t)  ”  0,  AQ  becomes  infinite  and  the  ray  method  fails.  In  the 
following  section,  we  shall  construct  uniform  asymptotic  expansions  to  remove 
these  anomalies. 


3.  Uniform  Asymptotic  Expansions 

The  main  idea  of  constructing  a  uniform  asymptotic  expansion  for  (4)  to 
(6)  is  based  upon  the  assumptions  that  a  reflected  ray  emerges  after  an 
incident  ray  reaches  an  anomaly,  and  the  behavior  of  the  solution  there  may  be 
asymptotically  approximated  by  the  solution  of  an  ordinary  differential 
equation  as  a  comparison  equation.  The  ansatz  is  the  following: 

♦  =  exp(-$0)[<J>mv(8U<;)  +  eU-V2>V'(eUC)]  ,  (15) 

Here  satisfies 

v“(8wo  +  -  82_2,jq2(c)v(bijc)  =  o  ,  ne) 


6 


{K*!1*)  -  K  76  .  [-♦! 1)  ve  +  ♦*2)Q2VO  +  KQ27?  •  [-V6*J2) 

Iz  Z  u  u  u 

-S1[-#J1,0t*#J2)«tfl2]  . 

<k*!2))  +  K[-^nvc  •  ve  +  ^2)Q2(7C)2J 

1  z  z  u  u 

-  K[-*‘2)<ve)2  +  tj^ve  .  v«  -  s1[-^2)et  +  ^nct]  , 


at  z  ■  0, 


at  z  ■  -h, 


+  S2["*O1>0t  +  *<2)ctQ21  =  0  ‘ 

*♦"’  *  s2‘-*52>9t  *  ♦o’’ct>  ■  0  • 


xf1 )  .  *(2>  _  0 
♦l«  *1z  0  • 


By  multiplying  (21)  by  ±Q  and  adding  the  resulting  equation  to  (20),  we 


obtain 


where 


«♦?,>,  -  <s/  -  Kl**)2!*?  • 


W*  =  -S*  -  -(6fc  +  QCt)  . 
K*  -  Vs*  -  V0  +  QVc  • 


Similarly,  (23)  to  (24)  imply 


«#*,♦  S2U±*J  *  0  at  e  =  0  * 


♦  ,  *  0  at  z  *  -h  • 

Iz 


Equations  (26)  and  (27)  are  of  the  same  form  as  given  in  the  ray  method 

expansion  [4].  Therefore,  it  follows  from  (8)  that 

+  ±2 
w1  *  D(kj  . 

The  equations  for  the  second  approximations  yield 
(X*2z )z  ~  -KCk*)2#*  +  «  (Kk)  +  2Kk  •  V*J  +  S<f%* 

+  S^Jt  -  D*S1#J2)Ct  -  K[2D±V8  .  VC^2)  -  D*( Vc)2^1’ 
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(29) 


-  Q<VC)V>J2)]  , 

**i.  ■  -v‘»?  -  ♦  s2D^o2,{t  «  *  -  °  • 

*2z  “  ^07h  *  ^  at  2  =  _h  '  (31) 

where 

D*  *  ±(Q'  +  QP)  •  (32) 


The  derivation  o£  (29)  is  a  little  tedious  but  straightforward  and  the  details 
are  omitted  here.  We  integrate  (29)  from  z  =  -h  to  z  =  0,  and  make  use  of 
(30),  (31)  to  obtain 

-S^^U^O.t)  -  S2*Jt  -  tjx^y.-hjk*  •  Vh  +  S2D±Ct«J2)  , 

=  /°hI(S10)±  -  Ktk*)2)**  +  •  (Kk*)  +  2Kk*  •  V(j,J  +  S1  ^gfc]  dz 

-/®h(D±S1Ct*J2)  +  KlD^Ve  .  VC^2)  -  D±(VC)a*J1>  -  D+Q(VC)2^J2>]}dz.  (33) 


Let 


*  R(C)A*,  2RR*  -  R2D+Q_1 


(34) 


It  is  shown  in  Appendix  A  that  (33)  may  be  reduced  to  the  same  transport 
equation  (13) 

[I(Aj)2)t  +  V  •  /°h  2K(aJ )2k±dz  +  2N±(A*)2  =  0  ,  (35) 

and  by  integration  along  a  ray,  we  have 

I(A*)2J*(t)  exp(-/*  2N±I*1dt)  =  constant  , 

or  by  (34), 


I(R)”2(^^)2J±(t)  exp ( -/5  2N±l“1dt)  =  constant  .  (36) 

0  *0 

Based  upon  (32),  (34)  and  (36),  we  may  choose  suitable  P,  Q  and  R  so  that 
tend  to  a  finite  limit  at  an  anomaly. 
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4 .  Applications 

First  we  consider  a  line  of  aero  depth  T  :  h  ■  0,  assumed  to  be 

J/> 

smooth.  We  choose  y  -  1,  Q(C)  -  1.  since  A  behaves  like  C  2  at  r,  as 

1/, 

shown  in  the  Appendix  B,  we  choose  R  -  v2  *  Then  by  (32),  (34), 

2  -1 

P  ( C )  [ log  (R  /Q)]'  ■  C  •  The  comparison  equation  (16)  becomes 

v(8c)  +  (ecrVieo  -  v<ee)  -  o  .  (37) 

The  solutions  of  (37)  may  be  expressed  in  terms  of  modified  Bessel  functions 
of  zeroth  order.  We  choose  V(0t)  *  IQ  ( 60  so  that  V( 80  is  regular  at 
C  -  0  which  corresponds  to  h  y  0.  In  Appendix  B,  it  will  be  shown  that 
indeed  tend  to  a  finite  limit  as  h  ♦  0.  From  (15)  to  (19), 

♦  -  exp(-86)t*(1)I0(8O  +  ♦(2)1'(80]  ,  (38) 

0  -  (S+  +  S_)/2,  t  -  <S_  -  S+)/2  , 

♦  (1>  *  if  ♦  D/2,  ♦(2>  -  ♦*  -  r)/2  . 

If  we  have  two  lines  of  zero  depth  in  a  fluid  region,  corresponding  to 

Vo 

C  •  0,  C  *  we  may  choose  V  *  0,  Q(0  *  1,  R(0  **  [CtC^-O]  2  .  Then  it 

follows  that  P(0  -  [(((,-<)]  ,/[C(C1-0] .  (16)  now  becomes 

[((C^OVMO)'  -  82C(C1-0V(0  -  0  . 

At  a  strictly  convex  caustic,  we  follow  Ludwig  [10]  and  choose  y  *  2/3, 

Vo  V, 

Q  -  O2  ,  p  »  0.  From  (32),  (34),  it  follows  that  R  “  C  4  and  (16)  becomes 

V"(82/3C)  -  02/3CV(82/3C)  -  0  , 

2/3  2/3 

which  is  the  Airy  equation.  We  may  choose  V(6  O  *  Ai(6  O  to  construct 
a  uniform  asymptotic  expansion. 

In  case  we  have  a  line  of  zero  depth  corresponding  to  C  ■  0  and  a 

i 

caustic  corresponding  to  C  ■  we  take  Q  ■  (Cj-O^  ,  P  »  C  ,  y  *  0  and 

Vo  V, 

R  *  ?/2(C1“?)/4«  Then  (16)  becomes 

V"<C)  +  c“V(C)  -  02(C1-C)V(C)  -  0  . 
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Finally,  if  there  are  two  caustics  corresponding  to  C  =  0  and  C  - 
(16)  assumes  the  form 

v-K)  -  e2c(?1-;)v  =  0  , 


•  1 ,  thi 


and  R  =  C/4  (5,-5  )1/4  . 


Appendix  A.  Derivation  of  the  Transport  Equation 
Here  we  show  that  under  the  conditions  (34),  (33)  is  reduced  to  (13).  By 
comparing  (33)  with  (12),  we  see  that  the  only  difference  between  them  is  the 
extra  term  on  the  right-hand  side  of  (33).  We  multiply  both  sides  of  (33)  by 
2$*,  and  make  use  of  (34)  and  the  derivation  of  (13)  from  (12)  to  obtain 

R2{[l(Aj)2]t  +  V  •  J®h  2K(aJ)2  k*dC  +  2N±(aJ)2} 

«  2RR,{-[I(Aj)23Ct  -  J°h  2K(aJ)2  i^dC  •  7? 

±  QaJ  !°_h  [SlCtAj2,?t  +  K(2V0  •  7«Vq2)  -  (Vc^A*1*  +  Q(7C)2Ag2>)]dC 

x  to) 

±  2QaJaJ  ,  (A. 1 ) 

where 

aJ1J  -  <A+  +  A‘)/(2Q),  A<2)  =  (Aj  -  A^)/2Q)  .  (A. 2) 

If  we  can  show  that  the  right-hand  side  of  (A.1)  vanishes,  then  (A.1)  is 
reduced  to  (13).  To  this  end,  we  first  take  the  upper  sign  in  the  right-hand 
side  of  (A.1),  and  by  (19)  we  replace  70  by  (k  +  +  k  )/2,  7j;  by 
(k  -  K  +)/(2Q)  and  by  (w+  -  w  )/(2Q).  Then  by  (A. 2)  we  collect  the 

coefficients  of  A*  and  aq,  which  vanish  identically  because  of  (28). 
Similarly,  the  right-hand  side  of  (A.1)  corresponding  to  the  lower  sign  also 
vanishes  by  the  same  derivation.  Thereore,  (35)  and  (36)  follows. 
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Appendix  B.  Regularity  at  a  line  of  zero  depth 

To  be  definite,  we  assume  that  K,  S1  and  S2  are  smooth,  bounded  and 

positive  at  F  j  h  =  0.  By  (8)  and  (9)  and  the  expression  for  N(x,y,t)  in 

(13),  it  is  easily  shown  that  N  =  0(h)  near  T.  since  I  >  0,  we  see  from 

(34)  and  (36)  that  behave  like  R[J*(t)]  ^ .  Since  R  =  ?/2  ,  we  show 

±  +  -3/, 

J  (t)  =  0(0  so  that  R[J  (t)]  2  tend  to  a  finite  limit  as  C  ♦  0.  For 

convenience,  a  coordinate  system  (t,£,n)  is  chosen  where  £,  n  are 

respectively  distances  along  T  and  in  the  normal  direction  of  T  toward  the 

fluid  region.  Since  £,  n  are  independent  of  t,  the  Jacobian  of 

transformation  from  t,  £,  n  -  system  to  t,  x,  y  -  system  is  given  by 

=  1  “  Kn  '  <B.1) 

where  k  is  the  curvature  of  T.  in  terms  of  t,  £  and  n,  (28)  becomes 

S?(1-Kn)”2  +  s2  *  -D~V  (B.2) 

£  n  t 

where  we  drop  the  ±  signs.  Based  uponn  the  assumptions  on  K,  S1  and 

S2,  it  follows  from  (9)  that 

D  =  /®  K  dC/I  =  DQh  +  0(h2 )  ,  (B.3) 

where  D0  ?  0  at  T.  We  further  assume  3h/3n  =  Y  >  0  at  T  and 
2 

h  =  Yn  +  0(n  ).  Therefore,  we  can  express  (A. 3)  as 

D  *  nf (£,n)  (B.4 ) 


where  f(£,0)  >  0.  From  (A. 2)  and  (A. 5)  we  have 

nS?(1-Kn)“2  +  (n/2S  )2  =  -ff~V  .  (B.5) 

s  n  t 

I/. 

Now  we  introduce  a  new  variable  h  =>  n2  ,  and  in  terms  of  it,  (A. 5)  becomes 


2  2  -2  2  -1 

n  s^O-ien)  +  s^/4  =  -f  st  , 


(B.6) 


sn  =  ±2(-f_1st  -  n2s|(i-icn2)_2]1/2  . 

(B.7 )  may  be  solved  by  the  method  of  characteristics  if  we  assume 


(B.7 ) 


uj  -  >  0  and  S  is  sufficiently  smooth,  on  h  *  0,  t  >  0.  It  is  not 


difficult  to  show  from  the  characteristic  equations  that  the  Jacobian  of 


transformation  from  the  ray  coordinates  t,  t0,  to  the  t,  C,  n  - 
coordinates  is  given  by 

J(r£t6,E  )  "  [“0f(€#°)]1/2  »  (B.8) 

where  a 1  =  tQ,  «2  *  is  a  point  on  r\  -  0,  t  >  0.  The  Jacobian  in  (36) 
may  be  expressed  as 


J*(t) 


(B.9) 


From  (A.1),  (A. 8)  and  (A.9),  we  see  that 

J*(t)  -  +2[o>0f(5,0)]1/2n  .  (B.  10 ) 

Now  from  (A. 7 ) , 

sn<t,C,n)  -  ±2[«0f(C,0)]1^  +  o(n2)  . 

By  integration, 

s*(t,C,n)  -  s0  +  2 («0f (5,0)}^  i)  +  o(n3)  , 

and  it  follows  from  (19)  that 

5  -  (s’  -  s+)/2  -  2[»0f(5,o)]1/2 n  +  o(n3)  .  (b.id 

By  (A. 10)  and  (A. 11),  we  have 

J*(t)  ~  +C  , 

and  tend  to  a  finite  limit  as  C  ♦  0. 
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